clear; clc; close all;

% 固定参数
a = 5; b = 1; P = 1; Q = 1;
x0 = 1; y0 = -1;
N_total = 1500;  % 迭代总步数
N_transient = 1000; % 丢弃暂态步数
delta0 = 1e-8;       % 初始扰动大小

% q范围和采样
q_list = linspace(0, 2, 200);
M_q = length(q_list);
MLE_vals = zeros(1, M_q);

for idx = 1:M_q
    q = q_list(idx);
    
    % 生成参考轨迹
    [x_ref, y_ref] = SCLMM(q, a, b, P, Q, x0, y0, N_total);
    
    % 初始化扰动向量（单位向量）
    v = delta0 * [1; 0];
    
    % 初始化扰动轨迹起点：参考轨迹初始点加扰动
    X_pert = [x_ref(1); y_ref(1)] + v;
    
    sum_ln = 0;      % 累计扰动扩展对数和
    count = 0;       % 有效重正化次数
    
    for n = 1:N_total-1
        % 参考点和扰动点
        X_ref = [x_ref(n); y_ref(n)];
        % 计算雅可比矩阵（用参考点）
        J = computeJacobian(X_ref, a, b, P, Q);
        
        % 利用雅可比矩阵演化扰动向量
        v = J * v;
        
        % 并行计算下一步扰动轨迹点
        X_pert = X_ref + v;
        
        % 每隔10步进行一次重正化（可调）
        if mod(n, 10) == 0
            norm_v = norm(v);
            sum_ln = sum_ln + log(norm_v / delta0);
            v = (delta0 / norm_v) * v;  % 重正化扰动向量
            count = count + 1;
        end
    end
    
    % 计算最大李雅普诺夫指数（平均扰动扩展率）
    MLE_vals(idx) = sum_ln / (count * 10); % 除以总迭代步长
    
    fprintf('q=%.3f, MLE=%.5f\n', q, MLE_vals(idx));
end

% 绘制MLE曲线
figure('Position',[100 100 900 600]);
plot(q_list, MLE_vals, '-r', 'LineWidth', 1.5);
xlabel('q', 'FontSize', 14);
ylabel('Maximum Lyapunov Exponent (MLE)', 'FontSize', 14);
title('Fig.4(b) 最大李雅普诺夫指数（Wolf方法）', 'FontSize', 16);
grid on;
box on;


% ======================
% SCLMM分数阶模型函数（你提供）
function [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N)
    x = zeros(1, N);
    y = zeros(1, N);
    x(1) = x0;
    y(1) = y0;

    w = zeros(1, N);
    for k = 1:N
        w(k) = exp(gammaln(k + q - 1) - gammaln(k));
    end

    delta_x = zeros(1, N);
    delta_y = zeros(1, N);

    for n = 2:N
        delta_x(n-1) = a * sin(P * y(n-1)) * sin(P * x(n-1)) - x(n-1);
        delta_y(n-1) = b * cos(Q * x(n-1)) * cos(Q * y(n-1)) - y(n-1);

        sum_x = 0;
        sum_y = 0;
        for i = 1:n-1
            sum_x = sum_x + w(n - i) * delta_x(i);
            sum_y = sum_y + w(n - i) * delta_y(i);
        end

        x(n) = x0 + sum_x / gamma(q);
        y(n) = y0 + sum_y / gamma(q);
    end
end


% ======================
% 雅可比矩阵计算函数（你提供）
function J = computeJacobian(X, a, b, P, Q)
    x = X(1);
    y = X(2);

    J = zeros(2,2);
    J(1,1) = -a * P * cos(P * x) * sin(P * y);
    J(1,2) = -a * P * sin(P * x) * cos(P * y);
    J(2,1) = b * Q * sin(Q * x) * cos(Q * y);
    J(2,2) = b * Q * cos(Q * x) * sin(Q * y);
end
